/*
Este programa implementa o algoritmo de Monte Carlo, como
apresentado em http://www.eveandersson.com/pi/monte-carlo-circle
Utiliza-se a biblioteca para aritmética de precisão GMP. A precisão mínima utilizada
em todas as variáveis relevantes do programa são tomadas como input. Neste
programa, a precisão é interpretada como o número de casas decimais que deve ser
impresso (não necessariamente corretos).
Realizam-se 10^precisão iterações.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <time.h>

int main(int argc, const char* argv[])
{
   //Declaração de variáveis para abertura de arquivo e medida de tempo
   FILE *fp;
   char filename[50];
   struct timeval ti, tf;
   long int sec, usec;

   //Início da contagem de tempo
   gettimeofday(&ti, NULL);

   //Declaração das variáveis utilizadas nos cálculos
   mpz_t iterations, m, n;
   mpf_t pi, aux;
   int precision;
   double x, y;

   //Inicialização de variáveis com a precisão especificada
   precision = atoi(argv[1]);
   mpf_set_default_prec(precision);
   mpz_init(iterations);
   mpz_init(m);
   mpz_init(n);
   mpf_init(pi);
   mpf_init(aux);
   mpz_ui_pow_ui(iterations, 10, precision + 1);

   //Loop principal
   while(mpz_cmp(n, iterations) != 0)
   {
      //Amostra um ponto no quadrado unitário
      x = ((double) rand() / (double) RAND_MAX);
      y = ((double) rand() / (double) RAND_MAX);
      
      //Incrementa o número de iterações
      mpz_add_ui(n, n, 1);
      
      //Se amostra estiver dentro do círculo unitário, incrementa m
      if(sqrt(x*x + y*y) < 1.0)
         mpz_add_ui(m, m, 1);
   }

   //Cálculo final de pi (pi =~ 4*m/n)
   mpf_set_z(pi, m);
   mpf_mul_ui(pi, pi, 4);
   mpf_set_z(aux, n);
   mpf_div(pi, pi, aux);
   

   //Término da contagem de tempo
   gettimeofday(&tf, NULL);
   sec = tf.tv_sec - ti.tv_sec;
   if(tf.tv_usec < ti.tv_usec)
      sec--;
   usec = (tf.tv_usec - ti.tv_usec + 1000000) % 1000000;
   printf("Tempo decorrido: %ld.%06ld s. Iterações: ", sec, usec);
   gmp_printf("%Zd\n", iterations);
   
   //Salva o resultado para arquivo
   sprintf(filename, "pi/montecarlo-seq-%d", precision);
   fp = fopen(filename, "w");
   gmp_fprintf(fp, "%.Ff\n", pi);
   fclose(fp);

   return 0;
}
